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MONOENERGETIC  NEUTRAL  PARTICLE  TRANSPORT 
IN  SEMI-INFINITE  MEDIA 


B.  Ganapolf 

Department  of  Nuclear  Engineering 
University  of  Arizona 


Abstract-  A  large  class  of  half-space  transport  problems  for  neutral  particles 
is  considered  including  single  and  adjacent  hcilf-space  geometries.  Localized 
and  distributed  sources  emitting  particles  either  in  a  particular  direction, 
isotropically  or  with  a  given  angular  distribution  are  considered.  All  so¬ 
lutions  are  based  on  the  solution  to  the  well-known  eilbedo  problem.  A 
numerical  Laplace  transform  inversion  is  used  to  generate  highly  accurate 
benchmark  solutions  which  can  be  used  to  assess  the  error  incurred  by  nu¬ 
merical  transport  algorithms. 


I.  INTRODUCTION 

Particle  transport  theory  had  its  inception  with  the  establishment  of  the 
Boltzmann  equation  by  Ludwig  Boltzmann  in  the  late  nineteenth  century^. 
In  its  most  general  form,  the  equation,  which  serves  as  the  basis  for  many 
fields  of  physics,  has  evaded  analytical  solution.  Only  for  special  cases,  i.e., 
the  linecir  or  linearized  forms  and  simplified  non- linear  collision  models,  have 
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analytical  solutions  been  obtained.  With  today’s  computational  power,  how¬ 
ever,  the  numerical  solution  to  the  Boltzmann  equation  for  physically  inter¬ 
esting  cases  in  the  linear  and  non-linear  formulations  can  be  obtained  rela¬ 
tively  easily.  Various  methods  such  as  discrete  ordinates,  Legendre  and  other 
polynomial  expansions,  the  Fn  method,  direct  simulation  Monte  Carlo  and 
other  approximate  methods  have  proven  to  provide  adequate  numerical  so¬ 
lutions  for  certain  applications.  These  methods,  however,  contain  numerical 
approximations  in  various  forms.  The  effect  of  these  numerical  errors  usually 
goes  unassessed  and  their  influence  on  the  desired  result  remains  uncertain. 

In  1953,  the  first  comprehensive  tool^  attempting  to  provide  standards 
to  allow  the  assessment  of  numerical  errors  associated  with  approximate  so¬ 
lutions  to  the  linear  neutral  particle  transport  equ^ttion  was  published.  This 
compilation  is  a  collection  of  analytical  solutions  complete  with  numerical  e- 
valuations  (some  of  which  are  unfortunately  inaccurate)  which  could  be  used 
as  standards  or  benchmarks  to  which  approximate  solutions  could  be  com¬ 
pared.  Infinite  medium  problems  were  emphasized  in  this  collection  because 
standard  numerical  methods  at  that  time  were  considered  reliable  enough  to 
adequately  evaluate  the  analytical  representations  obtained  for  this  class  of 
problems. 

It  is  now  the  nineties,  almost  50  years  since  the  first  analytical  bench¬ 
mark  compilation  appeared.  During  that  time  there  has  been  enormous 
progress  made  in  the  areas  required  for  the  generation  of  transport  bench¬ 
marks,  namely,  in  analytical  and  numerical  methods  development  and  com¬ 
putational  power.  This  report  is  concerned  with  the  development  of  the 
components  required  specifically  for  the  generation  of  transport  benchmarks 
in  one- dimensional  semi-infinite  (half-space)  plane  geometry. 

A  fundamental  problem  in  neutral  particle  transport  theory  concern- 
s  the  collision  of  neutral  particles  with  stationary  scattering  centers  and 
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the  subsequent  streaming  between  collisions  in  a  semi-infinite  medium.  The 
semi-infinite  homogeneous  medium  is  infinite  in  extent  bordering  either  on 
a  vacuum  or  another  half-space  of  a  different  material.  In  general,  upon 
collision  with  a  scattering  center  (usually  a  molecule  or  nucleus),  a  particle 
will  lose  energy  elastically  or  inelastically.  For  our  purposes,  we  will  consider 
only  monoenergetic  neutral  particle  transport  which  ignores  inelastic  scat¬ 
tering  and  assumes  heavy  scattering  centers  that  can  absorb  an  unlimited 


amount  of  momentum.  For  this  case,  neutral  particles  scatter  and  move  be¬ 
tween  collisions  with  a  constant  velocity.  Energy  transfer  considerations  are 
usually  handled  through  the  multigroup  formulation  of  which  monoenergetic 
transport  is  a  special  case.  In  addition  to  the  assumption  of  monoenergetic 
transport,  we  will  assume  that  the  neutral  particle  leaves  a  collision  with  no 
preferred  direction  which  is  the  assumption  of  isotropic  scattering. 

Half-space  transport  problems  first  appeared  in  connection  with  photon 
transport  in  planetary  atmospheres  in  the  form  of  the  Milne  problem^"®.  For 
this  case,  a  source  embedded  deep  within  the  interior  of  a  star  continuously 
emits  an  unlimited  supply  of  photons,  a  portion  of  which  exit  the  star’s 
atmosphere.  The  curvature  of  the  atmospheric  surface  can  be  neglected  in 
comparison  to  the  photon  mean  free  path  which  allows  a  plane  geometry 
approximation  to  be  used.  Variations  based  on  the  Milne  problem  soon 
appeared  in  the  literature  as  part  of  neutron  transport  theory®  where  the 
source  was  moved  to  the  free  surface  (called  the  albedo  problem).  Because 
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One  of  the  major  themes  of  this  investigation  is  the  establishment  of  half¬ 
space  solutions  based  solely  on  the  solution  to  the  single  half-space  albedo 
problem  for  an  entering  beam.  It  is  remarkable  that  almost  all  solutions  to 
half-space  transport  problems,  including  those  for  two  dissimilar  half-spaces, 
can  be  obtained  from  this  one  basic  solution.  Another  pervasive  theme  found 
here  is  the  ease  with  which  truly  accurate  benchmarks  can  be  generated.  The 
key  to  the  accuracy  of  the  numerical  evaluation  is  the  numerical  Laplace 
transform  inversion  which  has  recently  been  developed  to  solve  the  types  of 
transport  problems  to  be  considered. 

Single  half-space  problems  including  the  albedo,  Milne,  and  spatially- 
distributed  sources  will  be  considered  in  §  II.  The  preferred  method  of  so¬ 
lution  is  the  Laplace  transform  and  its  inversion.  The  motivation  for  this 
choice  is  the  desire  to  obtain  useable  numerical  results  rather  than  simply 
analytical  representations.  The  solution  representation  for  the  angular  and 
scalar  (angularly  integrated)  fluxes  can,  in  general,  be  expressed  as  contour 
integrals  which  can  further  be  manipulated  by  analytic  continuation  if  the 
singularities  in  the  complex  plane  are  known.  Finding  these  singularities 
usually  requires  a  significant  computational  effort  which  is  avoided  in  the  in¬ 
version  procedure  that  has  been  developed.  The  numerical  inversion  is  also 
applied  to  adjacent  half-space  problems  which  are  considered  in  §  III.  Final¬ 
ly,  in  §  IV  the  numerical  implementation  and  demonstration  of  the  solution 
technique  is  presented. 
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II.  SINGLE  HALF-SPACE  SOLUTIONS 


A  single  half- space  bordering  a  vacuum  will  now  be  considered.  The 
medium  is  assumed  to  be  homogeneous  and  to  contain  scattering  centers 
that  scatter  neutral  particles  isotropically.  Physically,  the  assumption  of 
isotropic  scattering  means  that  the  direction  of  the  "outgoing”  particle  after 
collision  is  randomly  oriented  and  has  no  correlation  with  the  direction  of  the 
"incoming”  particle  initiating  the  collision.  In  addition,  as  discussed  in  the 
introduction,  the  scattering  interation  does  not  change  the  particle  energy  (or 
velocity).  Various  source  configurations  will  be  considered  including  surface 
and  spatially  distributed  sources  and  an  embedded  source  infinitely  far  from 
the  surface. 

The  single  half-space  problem  is  one  of  the  fundamental  problems  of 
transport  theory.  It  not  only  is  of  pedagogical  interest  but  also  has  practical 
significance.  In  the  non-destructive  analysis  of  materials  with  high  mass 
numbers,  the  reflected  beam  resulting  from  an  impinging  neutral  beam  can 
be  used  to  infer  material  properties.  The  monoenergetic  single  half-space 
problem  provides  a  first  approximation  to  the  reflected  beam  and  can  help 
guide  the  solution  to  the  inverse  problem  from  which  the  desired  scattering 
characteristics  of  a  medium  can  be  inferred. 


II. 1  The  Albedo  Problem 

For  this  problem,  the  source  configuration  is  assumed  to  be  particles 
impinging  on  the  free  surface  of  a  homogeneous  semi-infinite  medium.  If  any 
absorption  exists  in  the  medium,  then  infinitely  far  from  the  surface  the  flux 
will  vanish.  However,  if  no  absorption  exists,  then  the  flux  will  remain  finite. 
Thus,  for  the  following  albedo  problems  considered  for  a  variety  of  source 
angular  distributions,  the  condition  at  infinity  will  be 
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lim  <  cxD  ,  0  <  c  <  1 

X— ►(» 

where  ^(x,/i)  is  the  directed  or  angular  flux,  c  is  the  number  of  secondary 
(scattered)  particles  resulting  from  a  collision  and  is  explicitly  given  by 

c  =  (S,  +  . 

Sr  represents  the  cross  section  (fractional  probability  per  pathlength  of  trav¬ 
el)  for  process  r  where 

r  =  s  scattering 
r  =  a  absorption 
r  =  f  fission 

r  =  t  total  (St  =  Sa  +  Ss)  . 

Fission  is  included  for  the  case  of  neutrons  with  u  representing  the  number 
of  neutrons  produced  in  a  fission. 

The  tremsport  equation  to  be  solved  for  the  classic  albedo  problem  is 

'  d  c 

^(0,/x)  =  f{f^)  ,  M  >  0  (16) 

lim  ^(a;,/i)  <  oo  (Ic) 
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where  x  is  the  position  measured  in  units  of  r:ean  free  paths  (x  =  zHt,  z  = 
actual  position),  /i  is  the  particle  direction  cosine  (referred  to  as  the  direc¬ 
tion)  and  /(/i)  is  a  general  source  distribution  to  be  specified  in  the  following 
subsections. 


l.a.  Entering  Beam  Flux 

The  solution  to  the  albedo  problem  for  a  beam  source  serv'es  as  the  basis 
for  the  solution  to  problems  with  the  source  configurations  to  follow.  For  this 
reason,  the  solution  to  the  beam  albedo  problem  will  be  presented  in  detail. 
For  a  beam  source. 


/(//)  =  fio)  ,  (2) 

where  the  source  is  specified  to  be  in  direction  /io  (which  is  actually  a  "cone” 
of  incoming  particles  in  the  azimuthally  symmetric  plane  geometry  assumed). 


i.  Exiting  flux 

The  exiting  or  reflected  flux  is  likely  to  be  the  most  important  quantity 
to  be  determined.  It  also  serves  2is  the  nucleus  for  the  determination  of  the 
interior  scalar  flux  as  well  as  the  angular  flux.  A  sketch  of  the  determination 
of  an  analytical  expression  for  the  exiting  flux  will  now  be  given. 

The  procedure  begins  with  the  reformulation  of  eqs.(l)  by  following  the 
particle  along  its  trajectory  to  give  for  fj,  >  0 

7p{x,fj,;fio)  =  6{fi  -  ^  (3a) 

2^  Jq 


and 
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J  X 


where  the  scalar  flux  is  defined  as 


V>(x;/io)  =  J  d/iXx,/i';/io)  (4) 

and  the  dependence  on  fio  has  been  included  for  clarity.  Then  integrating 
eqs.(3a)  and  (3b)  over  on  [0, 1]  gives  the  following  integral  equation; 

^  /  di'£i(l  X  -  i'  |)r/;(x':^o)  (5) 

-  Jq 

where  Ei  is  the  exponential  integral  defined  by 


Ei{x) 


-L 


This  integral  equation  can  be  put  in  the  more  compact  notation 


(1  -  Ii)^(x';/io)  =  e 


where  the  operator  Li  is 


I,  =  r  dx'Er{\x-x'  |)(.)  . 
Jo 


Differentiating  eq.(6a)  with  respect  to  x  yields 


and  also  memipulating  eq.(6a)  gives 
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(8) 


£;,(x)  =  (l-I,)  /  ^V-(x';/)  . 

Jii 

From  eqs.(7)  and  (8),  it  is  seen  that  the  following  expression  is  in  the  null 
space  of  the  operator  1  —  Lx: 

4-  +  —  =  lv>(0;^o)  /  •  (9) 

_dx  /xoj  2  Jo  fJ' 

The  reciprocity  relation, 

fiip(0,-fi-,fio)  =  fiorpiO,-Ho;  ^J.)  ,  (10) 

obtained  from  the  self-adjointness  of  Lx,  can  be  introduced  into  the  result  of 
multiplying  eq.(9)  by  and  integrating  over  x  to  give  [after  liberal  use 

of  eq.(3b)  with  x  =  0] 

V’(o,-^;/io)  =  V^(0;/^o)V’(0;m)  ■  (H) 

2  ft  +  fio 

When  this  expression  is  introduced  into  the  following  representation  of  the 
scalar  flux  at  x  =  0  [eq.(4)  with  x  =  0] 

^(0;/xo)  =  1  +  /  (12) 

Jo 

and  defining  the  if— function  by 

if(/x)  =  V’(0;/i)  ,  (13) 


there  results  for  the  exiting  flux 
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2  /i  +  /io 


The  Chandresekhax  H — function"*  can  be  shown  to  satisfy  the  following  non¬ 
linear  integral  equation 


HM  =  1  +  |/<»(»‘)  / 


n  +  n' 


which  will  serve  as  the  basis  for  the  numericcil  evaluation  to  be  presented  in 

§IV. 

ii.  Scalar  flux 

An  expression  is  most  eeisily  obtained  for  the  interior  scalar  flux 
using  eq.(3b)  with  x  =  0 

V>(0,-/z;/io)  =  77-  /  dx'  .  (16) 

2/^  Jo 

Noting  that  if  //  is  extended  to  the  complex  variable  1/s,  eq.(16)  becomes  a 
Laplace  transform  with  the  inversion 


ip{x;no)  -  [V-CO, -l/s;/io)/s] 

c 


where  *  is  the  inverse  Laplace  transform  operator 


■1  f-r-^ioo 


Since  an  analytical  expression  for  the  exiting  flux  already  exists  [eq.(14)], 
eq.(17)  becomes 


misY 

1  +  HqS 


(18a) 


with 


^(1/5) 


1 


<i/i 


>  HW  • 

1+^'s 


(186) 


as  obtained  from  eq.(15).  Eq.(18a)  can  further  be  expzmded  by  analytical 
continuation  into  the  complex  s-plane  leading  to  contributions  from  pole 
and  branch  point  singularities  if  desired.  This  will  not  be  done  here  since 
the  numerical  inversion  given  in  §  IV  will  be  applied  directly  to  eq.(18a). 


iii.  Angular  flux 

The  most  straightforward  way  to  obtain  the  angular  flux  is  to  apply  the 
Laplace  transform  in  x  to  the  original  transport  equation  [eqs.(l)]  to  give 

—  c  — 

{1  +  ■  (19) 

From  eq.(17),  by  application  of  the  Laplace  transform,  we  have 

^(5;a^o)  =  -V’(0,-l/5;/io)/5  ,  (20) 

c 

resulting  in  the  transformed  angular  flux  from  eq.(19): 


rp{s,fi-,no)  =  — — [x(j{0,-l/s-,no)/s  +  /iV>(0,//;^o)]  • 

I  +  IJiS 


Thus,  for  fj,  >  0 
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[(1  +  fj.s){l  +  f^os) 


(22a) 


and 


1 

Hifi)  n 

1  -  fis  \  1  +  Hqs 

1  +no/nj_ 

(226) 


The  angular  flux  is  therefore  given  by  the  inversion 

V>(a:,/z;/io)  =  [^(5,  ^o)] 


(23) 


Again  further  manipulation  is  possible  using  analytical  continuation  into  the 
complex  s-plane  but  is  not  required  in  order  to  obtain  reliable  numerical 
results  as  will  be  shown. 


l.b  Entering  Isotropic  Flux 

For  an  isotropic  source,  the  boundary  condition  is 

V’(0,/i)  =  1  ,  /i  >  0  .  (24) 

where  the  entering  flux  has  been  normalized  to  1. 


i.  Exiting  flux 

The  exiting  flux  for  this  case  is  obtained  by  integrating  the  result  for 
the  albedo  problem  over  fiQ  on  [0, 1] 
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(25) 


fi;  Ho)  , 


to  give 


iP{0,-h)  =  1  -  \/l  -  c 


(26) 


where  the  relation  given  by  eq.(A.2)  in  Appendix  A  has  been  used. 


ii.  Scalar  flux 

At  a;  =  0,  the  scalar  flux  is 

/■i  2 

V>(0)  =  /  diJioH(Ho)  =  -(1  -  \/l  -  c  )  (27a) 

Jo  c 

using  eq.(A.l).  For  a;  >  0,  the  scalar  flux  is  obtained  by  integrating  eq.(18a) 
over  Ho  and  using  eq.(A.2)  with  h  replaced  by  1/s  to  give 

[(1-VT^/7(1/s))/s]  .  (276) 


iii.  Angular  flux 

Integrating  eq,(21)  over  ho  gives 

^(«>/^)  =  — lV’(0,-l/s)/5  +  /^0(O,^)]  ;  (28) 

1  +  /iS 

and,  therefore,  for  /i  >  0  using  eq.(26)  as  is  and  also  with  h  replaced  by  1/s, 
we  have  the  transforms 
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-cHil/s)]  (29a) 


TT^s  +  ^ 

\  ~  {(1/^  -  m)  -  >/l  -  c  [H  (l/s)  /s  -  nH{^i)] }  (296) 


which  axe  to  be  inverted. 


l.c  Entering  Angularly  Distributed  Flux 

For  a  general  angularly  distributed  entering  flux,  the  boundary  condition 


V>(0,^)  =  /(/x)  ,  for  /i  >  0  . 

This  case,  as  for  the  previous  cjise,  is  handled  by  noting  that  for  an  entering 
beam  the  flux  is  actually  the  ’’angular”  (partial)  Green’s  function;  thus, 

=  I  dfi’ .  (30) 

Jo 


i.  Exiting  flux 

Using  eq.(14),  the  exiting  flux  simply  becomes 

m  -m)  =  f  d^'-^fiix'W)  .  (31) 

Jo  M  r  /X 

ii.  Scalar  flux 

For  X  =  0, 
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(32a) 


JQ 


and  for  X  >  0  from  eqs.(30)  and  (18a) 


V>(x)  = 


(326) 


iii.  Angular  flux 

Integrating  eqs.(22)  over  /xq  on  [0, 1]  weighted  by  /  yields  for  /li  >  0 


V>(5,/X)  = 


1  fts 


fW 


c  H(l/s) 

I 


2(1  +  fis) 


i: 


dfi' 


1  +  fi's 


fifi'Wti')  (33a) 


2  1  -  fxs  Jo 

with  the  inversion  providing  the  angular  flux  as 

=  £"^[^(s,/i)]  . 


1  +  fi's 


(336) 


(34) 


II. 2  The  Milne  Problem 

The  Milne  problem  is  one  of  the  first  half-space  problems  ever  attempted. 
This  problem  found  application  in  the  determination  of  the  photon  current 
exiting  the  atmosphere  of  a  star  where  the  source  is  deeply  embedded  within 
the  star’s  interior.  For  this  case,  the  somrce  is  at  infinity  and  must  be  of 
infinite  strength  itself  in  order  to  support  a  current  infinitely  far  away  exiting 
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the  surface.  The  traditionad  condition  taken  at  infinity  is  that  the  flux  varies 
as  the  discrete  eigenmode  of  the  transport  operator 


</>(!, ^)  — ♦  ,  as  X  — >  oo  (35a) 

where 

=  I  ,  (356) 

2  /i  +  t'o 

and  vq  is  the  solution  to  the  dispersion  relation 

'  ■  I""  ■"  (^)  =  “  ■ 

The  transport  equation  to  be  solved  is  therefore 

^  ^  “  ly  dn'<t>{x,n')  (36a) 

=  0  ,  ^  >  0  (366) 

<^(x,/x)  — +  <^„o(— ^)e*/‘'®  ,  as  X  — »  oo  (36c) 

which  is  recast  into  a  more  convenient  form  by  the  substitution 

4>{x,n)  =  0^0 ( -/i)  e*/*'®  -  xl;{x,fi)  (37) 

yielding 
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(38a) 


d  c  f  I 

^'dx  ^  =2/1 

t/>(0,/i)  =  ,  fi  >  0  (386) 

lim  ip(x,fi)  =  0  .  (38c) 

r— *00 

It  is  this  half-space  problem  which  is  to  be  solved  based  on  the  solution  for 
the  general  angularly  distributed  entering  flux  found  in  §  II.l.c. 

2.a.  Exiting  flux 

Since  eqs.(38)  represent  a  half-space  problem  with 

fit^)  =  M-t^)  ^ 

eq.(31)  can  be  used  directly  to  give 

Employing  the  identity 

_ // _  _  1  y.  _  vq 

{y  -t-  y'){yQ  +  y')  y-  uq  [y’  +  y  y'  +  uo 

and  the  integral  equation  for  the  H—  function  [eq.(15)]  gives  the  following 
(after  some  algebra): 

^(0,-.)  =  KW  -  (40) 
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which  when  substituted  into  eq.(37)  (with  a:  =  0)  yields  the  desired  exiting 
flux 


m-f^)  =  o 


c  1/0  ^(m) 


2u(j-  n  H{uq) 


(41) 


2.b.  Scalar  flux 

From  integral  transport  theory,  as  found  in  §  II.  1, 


>A(0,  -/<)  =  ^ 

Jo 


(42) 


from  which  results  (as  previously  noted) 


^(x)  =  ■ 

c 


Thus,  form  eq.(40) 


'4’{x)  =  e 


—  _ 


HivoV"  [suo-l\ 


^0  r-l 


(43) 


and  from  eq,(37),  integrated  over  /r  noting  that 


dfi'M-fi')  =  1 , 


there  results 


=  litf" 


g(iA) 

Sl/Q  —  1 


(44) 


2.C.  Angular  flux 
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As  before,  the  spatiaJ  transform  of  0(x,/i)  is  given  by  eq.(28)  with  xp 
replaced  by  (f>  leading  to  for  //  >  0 


7,  X  ^  £ _ £0 _ Hjl/s) 

2(l  +  /is)(si/o-l)  H{uo) 

It  \  ^ 

2  H{vq){1  -  tis) 

with  the  angular  flux  again  given  by 

=  C~^  [<^(5,/i)]  .  (46) 


Hjlh)  _ 

SVq  —  1  Uq  —  fi\ 


(456) 


II.3  Spatially  Distributed  Sources 

The  half-space  problem  with  a  spatially  distributed  isotropicedly  emit¬ 
ting  source  will  be  considered  in  this  section.  The  appropriate  transport 
equation  to  be  solved  is; 

=  l/  dn'<l>{x,fi')  +  (47a) 

=  0,  ^  >  0  (476) 

lim  (f>{x,fi)  <  oo  .  (47c) 


The  solution  will  again  be  obtained  using  the  results  of  §  II.  1. a. 

Starting  with  eqs.(l)  for  an  entering  beam  flux  and  separating  the  un¬ 
collided  from  the  collided  contribution,  we  have 

xI>{x,h\hq)  =  V’o(x,/i;/io)  +  xpc{x,fi;fio)  ;  (48) 
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Jind,  therefore, 


51  c  c 

1  ^c{x,h\Hq)  =  2  J  (49a) 


0c(x,/i;/io)  =  0  ,  /i  >  0 


(496) 


lim  Ido)  <  oo 


(49c) 


If  /io  is  extended  to  the  complex  p-plane  by  replacing  /io  by  -1/p  and  eqs.(49) 
are  operated  on  by 


1  f'l+ioo 

—  y  dp  Q(p)(-)  ,  (50) 

eqs.(49)  become  eqs.(47)  with  (from  the  uniqueness  of  the  solution) 

j  1  ry+ioo 

(t>ix,p)  =  — -  /  dp  Q{p)  \pc{x,p-,-l/p)  (51) 

ZTTt  C 

where  the  Laplace  transform  of  the  source,  Q{p),  is  assumed  to  only  have 
singularities  in  the  left  half-plane  for  simplicity.  It  is  this  relation  to  the 
beam  albedo  problem  that  will  be  exploited  in  order  to  derive  the  results  to 
follow. 


3.a.  Exiting  flux 

From  eq.(51)  with  x  =  0  and  p  replaced  by  — p  and  from  eq.(14)  with 
po  replaced  by  —1/p  (note  that  the  uncollided  contribution  is  zero  at  x  =  0), 


,  cH(p)H(-l/p) 

V»c(0,-p;-l/p)  = 

2  1  —  pp 


there  results 


m-f^)  = 


1  H(fi) 
27ri  2 


dp  Q{p) 


I-  pp 


For  the  special  case  of  an  exponential  source,  we  find 


Q(x)  =  c-“*  ,  a  >  0  , 


(54a) 


which  implies 


Qip)  = 


p  +  a 


(54b) 


Eq.(53)  becomes  upon  evalution  of  the  contour  integral  lff(—l/p)  has  singu¬ 
larities  in  the  right-half  plane] 

I'+a—- 


3.b.  Scalar  flux 

From  integral  transport  theory  applied  to  eqs.(47),  we  find 


dx'  +  ^  /  dx'  e-’’'l'‘Q(x') 

Jo  2^  Jo 


Thus, 


c  c 
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with  <f>(0,—l/s)  given  by  eq.(55);  for  the  particular  case  of  an  exponential 
source 


<f>{x) 


(58) 


3.C.  Angular  flux 

Taking  a  Laplace  transform  of  eqs.(47)  yields 


(f>(s,n) 


2  1  +  /is 


1  Qjs)  ^ 

21  +  fis 


1+  fis 


HO, It)  ■ 


(59) 


Since  from  eq.(57) 


-(^(0,-l/s)/s 

c 


$(£) 

c 


eq.(59)  becomes 


4>{s,n)  =  — [<(>{Q,-l/s)/s  +  ;  (60) 

1  +  fiS 

and,  therefore,  for  //  >  0 


<f>{s,n)  =  (f>{0,-l/s)/{s{l  +  ns))  (61a) 


<l>{s,-n)  =  — [<^(0, -l/5)/5  -  /i<^(0,-/x)]  .  (616) 

1  —  //s 


For  em  exponential  source  again  using  eq.(55)  for  /i  >  0 


(i>{s,n) 


1  Hil/a)H(l/s) 

2  (1  +  ns){s  +  a) 


(62a) 
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21^  TT^  "  ^TiT^ 


(626) 


with  the  angular  flux  recovered  from  eq.(46). 
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III.  ADJACENT  HALF-SPACE  PROBLEMS 

For  some  physical  situations,  especially  in  electron  transport,  interfa¬ 
cial  effects  between  two  media  are  important.  For  this  reason,  the  adjacent 
half-space  configui.^tion  will  be  investigated.  The  geometry  consists  of  two 
adjacent  half-spaces  of  differing  scattering  and  absorbing  properties.  Severzil 
of  the  source  configuratons  of  the  previous  section  will  be  considered.  .4gain 
the  2malytical  solution  are  based  on  the  solution  to  the  beam  albedo  problem 
of  §  II.  1. a. 


III.l  Interfacial  Flux  Source 

The  spatial  domain  is  extended  to  negative  positions  measured  in  units 
of  the  mean  free  paths  of  each  medium.  The  transport  equation  to  be  solved 
in  each  medium  (j=l,2)  for  an  interfacial  flux  source  is: 


'  9  1 

+  'J 

c  , 

<t>j(x,n)  =  y  y  dfj,'(f>j{x,n') 

(63a) 

+  /(/i) 

(636) 

lim  0i(x,/i)  <  oo  ,  lim  <  c»  . 

(63c) 

Z—*00  X— ♦  — OO 


The  half-spaces  eire  positioned  as  shown  in  Fig.l.  Several  flux  source  condi¬ 
tions  [/(/i)]  at  X  =  0  will  be  considered  in  the  following  sections. 


Medium  2 


Medium  1 


Figure  1 


24 


l.a  Beam  Interfacial  Source 

The  analysis  follows  for  a  beam  flux  directed  into  medium  1: 

fin)  =  6(ti  -  /io)  ,  /io  >  0  .  (64) 


i.  Exiting  flux 

The  analysis  begins  by  reformulating  eqs.(63)  in  terms  of  the  collided 
and  uncollided  flux  components.  Thus,  if 


then 


<^j,o(x,M-,Mo)  =  po)  (66) 


and 


+ 


1 


02,c(O,//;/io) 
lim  (^i,c(a:,/z;^o)  <  oo 


c 

-  ~2  J  ^  <^/^V>,c(x,/i';//o)  + 


=  ^!>i,c(0,/i;/io) 


,  lim  <f>2,c{x,fi-,no)  <  oo  . 


(67a) 

(676) 

(67c) 


The  results  of  Prinja^  axe  central  to  the  analysis  to  follow.  The  solution- 
s  found  in  ref.  7  are  beised  on  the  solution  to  an  integral  equation  formed 
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by  considering  each  half-space  separately  with  eq.(63b)  serving  as  the  con¬ 
nection.  Eq.(30),  with  an  appropriately  defined  /(/i)  for  each  half-space, 
then  provides  a  set  of  integral  equation  to  be  solved.  The  solution  to  these 
equations  [letting  p  =  l/po  in  eq.(48)  of  ref.7]  for  the  exiting  flux  is 


<i>2iO,p;pQ) 


C1C2 

2(c2  -  Cl) 


Cl  Hi{po)Hi{p)'  1 

C2  H2{ij.q)H2{p)\  I  +  pIpq 


(68a) 


C1C2 

2(c2  -  Cl) 


1 


Hiipo)H2{p))  1 

H2{po)Hi{p)\  1-  fi/fio 


(686) 


In  addition,  a  set  of  integral  equations  for  p.\  po)  cam  be  directly 
established  by  considering  the  incoming  flux  for  each  half-space  to  be  known. 
For  po  >  0  and  using  eq.(31)  in  medium  1  with 


f{p)  =  8{p-pq)  -f  <i>2{0,p-,PQ)  (69) 


amd  in  medium  2  with 


/(/^)  =  <^i(0,-^;^o) 


(70) 


gives 


2  p-\-  pq 

+  f'  (71o) 

*  Jo  y- 

f  (716) 

^  Jo  y ''  y 
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This  set  of  integral  equations  is  numerically  convenient  as  an  alternative  so¬ 
lution  especially  for  n  =  fiQ  and  ci  —  C2  when  eqs.(68)  become  indeterminant 
and  therefore  difficult  to  evaluate.  Also  eqs.(71)  can  be  used  to  generate  the 
interfaciaJ  distributions  for  a  general  interfacial  source  with  ease: 


WO,-*-)  =  f  V-^»2(/-')<A.(0,-/)  . 

2  Jo  11  +  fi' 


(72a) 


(726) 


ii.  Scalar  flux 

When  eqs.(63)  for  the  beam  flux  source  are  reformulated  using  integral 
transport  theory,  there  results  at  x  =  0 


poo 

</»i(0, -/x;/io)  =  /  dx'  e~^ 

Jo 

pOO 

<?^2(0,/i;/io)  =  ^  dx'  e~^  ^^^2{-x']fio) 
2^^  Jo 


which  implies 


(73a) 


(736) 


2 

<^i(x;/xo)  =  — £JM^i(0,-1/s;/^o)/5] 
Cl 

2 

hi-x;  Ho)  =  —C'^[(f>2{0,l/s]fio)/s] 

C2 


(74a) 


(746) 


where  (f>j{0,+l/s]fio)  are  given  by  eqs.(68)  or  (71). 


iii.  Angular  flux 

From  the  Laplace  trauisform  of  eqs.(63)  for  0  <  i  when  j  =  1  and  x  <  0 
when  j  =  2  respectively,  we  have 

—  1  r  c '  — 

+  /i<^j(0,/i;/io)J  (75) 

where 

Ho) 

with 

2 

=  — <^j(0,t1/£;^o)/5  (76) 

where  =F(— ,  +)  is  for  j  =  1,2  respectively,  and  eq.(75)  becomes 

^j{s,H]Ho)  =  yT — [?^j(0>Tl/«;^o)/s  +  /i(A>(0,/i;^o)]  •  (77) 

1  ”t“  fJiS 

The  inversion  is  then  given  by  eq.(46)  to  provide  the  angular  flux  in  each 
half- space. 

III. 2  Milne  Problem 

Because  the  adjacent  half-space  Milne  problem  is  well  documented  in 
the  literature®,  only  the  results  will  be  presented  in  this  section. 


=  /  dx  e  "0i(x,/i;|xo) 

Jo 

_  f°° 

=  dx  e  *^<f>2{-x,H;Ho) 
Jo 
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I 

i.  Interfacial  flux 

For  this  case  the  source  is  assumed  to  be  at  infinity  in  medium  1  yielding 


<^i(0, -/x) 


Cl  ^'0,1  Hiifi)  H2(uo,i) 

2  uo,i  -  n  Hi{y.)  Hx{vo,i) 


<^2(0, /x) 


( 


C2  t^O.l  H2{vq,\) 

2  t'o.i  +  Hi(fi)  Hi{uo,\) 


where  1/0,1  satisfies  eq.(35c)  with  c  replaced  by  ci. 


(78a) 

(786) 


ii.  Scalar  flux 

Again  form  integral  transport  theory,  one  finds 

4>i{x)  =  — [<^i(0,-l/s)/s]  (79a) 

Cl 

<t>2{-x)  =  — [<?i2(0,  l/s)/s]  .  (796) 

C2 

with  ^j(0,  q^l/s)  given  by  eqs.(78). 

iii.  Angular  flux 

As  in  the  previous  case,  eqs.(77)  apply  for  ;tx  >  0 

^  [«Aj(0,TlA)A  +  Ai<^>(0,/x)]  ,  (80) 

i  “T  jJtS 

and  the  angular  flux  is  obtained  by  Laplace  transform  inversion. 
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II.3  Spatially  Distributed  Source 

Here,  we  assume  a  spatially  distributed  isotropic  source  exists  in  the 
right  half-space  (j  =  1)  only.  Thus,  the  transport  equation  to  be  solved  is 


<i)j{x,^l)  =  y  dn'<j>j{x,n')  -I-  (81a) 

=  h{0,n)  (816) 

lim  <^i(x,/z)  <  oo  ,  lim  <i)i{x,n)  <  oo  .  (81c) 


For  hq  >  0,  the  transport  equations  for  a  beam  flux  source  for  the  collided 
contribution  are  given  by  eqs.(67).  As  for  the  single  half-space  case,  if 
is  replaced  by  —1/p  and  eqs.(67)  are  operated  upon  by  eq.(50)  where  Q{p) 
has  singularities  only  in  the  left-half  of  the  complex  p-plane,  then  it  can  be 
shown  that  for  j  =  1,2 

4>j{x,p)  =  /  dp  Q(p)<^j,c(x,p,-l/p)  (82) 

2m  Cl 

satisfies  eqs.(81).  The  exiting  scalar  and  angular  flux  for  this  case  will  be 
obtained  for  an  exponential  source  only  [see  eqs.(54)]. 


i.  Exiting  flux 

Since 


Q{p)  = 


1 

p  +  a 


the  contour  integration  in  eqs.(82)  yields  for  p  >  0 


(83) 
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(84a) 


4>2i0,ti)  =  (^2,c(0,/i;  l/a)/ci  (846) 

where  <^j,c(0,  1  /a)  for  j  =  1,2  are  given  by  eqs.(68)  or  (71)  with  replaced 

by  1/a. 


ii.  Scalar  flux 
As  in  §  Ill.l.a.ii  and  II. 3. b 


M^)  =  (85a) 

Cl  (Ci  ) 

(j)2{-x)  =  [<^2,c(0,1/5;  l/Qf)/5]  .  (856) 


iii.  Angular  flux 

The  transform  of  the  angular  flux  is  given  by 


1  +  ^5 


[(l)j,c{0,Tl/s;l/a)/s  +  /i(^_,,c(0,/i;l/a)]  ;  (86) 


and  finally  the  inversion  is  obtained  from  eq.(46)  in  each  half-space. 


31 


IV.  Numerical  Implementation  and  Demonstration 

At  the  heart  of  the  numerical  evaluation  of  the  exiting,  scalar  and  angu¬ 
lar  flux  lies  the  evaluation  of  the  H —function.  Currently,  there  exists  several 
methods  for  evaluating  the  //—function  including  the  evaluation  of 

-  the  analytical  representation  as  provided  by  the  Wiener-Hopf 
solution 

to  eq.(15) 

-  the  multiple  collision  solution  to  eq.(15) 

-  approximate  representions 

-  eq.(15)  directly  by  iteration. 

It  is  the  last  method  which  has  proven  to  be  the  most  reliable  as  well  as  most 
efficient. 

The  method  begins  with  the  discretization  of  the  fi  variable  using  the 
Gauss-quadrature®.  By  approximating  the  integral  in  eq.(15)  with  a  shifted 
Legendre- Gauss  quadrature  of  order  Lm  and  evaluating  the  resultant  equa¬ 
tion  at  the  quadrature  abcissae,  the  following  iteration  scheme  results: 


- 


Ht 


1  -1 


m=l 


(87) 


where 


if*  = 


(88) 


The  full  iteration  step  is  accomplished  via  satisfying  the  zeroth  moment: 


//*+*  = 


Oo 

do 


(89) 


with 
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(90) 


i'm  j 

m—l 

The  second  numerical  procedure  of  importance  is  the  numerical  Laplace 
transform  inversion.  The  numerical  inversion  is  obtained  through  the  follow¬ 
ing  steps: 

-  transformation  of  the  Bromwich  contour  formulation 

1  /•'ir+too 
A-.OO 

where  7  is  greater  than  the  largest  real  part  of  emy  singularity  of  the 
image  function  f(s),  into  the  cosine  integral 

fi^)  =  — T  /  dw  ^{f{'y  +  iw/t)]  cos{w)  (92) 

TTl  Jq 

-  reformulation  of  the  infinite  integral  into  an  infinite  series 

f{t)  =  ^(  —  1)*  /  dw  +  {w  +  ikn)/t)}  cos{w)  (93) 

k=o 

-  evaluation  of  the  integral  in  eq.(93)  with  the  Romberg  integration 
rule^° 

-  acceleration  of  the  convergence  of  the  infinite  series  through  either  the 
epsilon  or  the  Euler-Knopp®  acceleration  algorithm. 

Finally,  in  some  ceises,  it  is  desirable  to  solve  the  integral  equations 
given  by  eqs.(71)  rather  than  use  the  analytical  expressions  of  eqs.(68).  This 
is  simply  done  by  discretization  and  subsequent  matrix  inversion. 

Figures  2-6  axe  presented  as  a  demonstration  of  the  evaluation  of  the 
scalar  flux  for  a  selected  set  of  problems.  While  only  curves  axe  presented 
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here,  the  scalar  fluxes  for  all  cases  considered  in  this  report  can  be  obtained  to 
5-place  accuracy  or  better.  Thus,  the  numerical  Laplace  transform  inversion 
represents  a  true  analytical  benchmark. 

Figure  2  gives  the  flux  for  the  albedo  problem  for  a  normally  inci¬ 
dent  beam  flux.  The  number  of  secondaries  c  has  been  specified  to  be 
c  =  0.1(0. 2)0. 9.  Note  the  exponential  decrease  far  form  the  surface,  eis  one 
would  expect,  and  the  change  of  the  sign  of  the  derivative  with  c  at  x  =  0. 
Figure  3  provides  the  flux  for  the  same  variation  of  c  for  the  Milne  problem. 
This  time  the  flux  increases  exponentially  at  large  distances  and  deviates 
from  exponential  behavior  near  the  medium  surface.  Figure  4  shows  the  be¬ 
havior  of  the  scalar  flux  for  a  spatially  distributed  exponential  source  given 
by  eq.(54a)  with  a  =  0.0(0. 5)2.0.  For  a  =  0,  the  source  is  uniform  and  the 
flux  can  be  shown  to  asymptotically  approach  1/(1  —  c)  which,  in  this  case, 
is  10. 

Finally,  Figs. 5  and  6  give  the  scalar  flux  in  two  adjacent  half-spaces.  Fig¬ 
ure  5  depicts  the  flux  for  the  Milne  problem  for  a  variation  of  C2  =  0.1(0. 2)0. 9 
in  medium  2  and  c\  =  0.9  in  medium  1.  Note  the  strong  purturbation  of 
the  flux  near  the  interface.  Figure  6  gives  the  scalar  flux  for  a  uniformly 
distributed  source  (a  =  0  in  medium  1)  for  a  variation  of  C2  =  0.1(0. 2)0. 9. 
Again  the  scalar  flux  approaches  1/(1  —  cj)  asymptotically. 
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Fig. 2  Variation  of  c  for  albedo  problem  with  beam  source 


Fig.3  Variation  of  c  for  Milne  problem 
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APPENDIX  A 

Some  Useful  Relations  for  the  If-Function 

There  exists  several  useful  relations  satisfied  by  the  ff— function  which 
are  required  to  perform  /f— function  calculus.  These  include: 

A.l  Zeroth  moment 

1  o  _ _ 

=  -  [1  -  ]  (A.l) 

A.  2  An  integral  relation 

/i  =  f  =  1  -  VT^c  Hit,)  (A.2) 

APPENDIX  B 

Description  of  SEMIl  Program  for  the  Determination 

of 

Fluxes  in  a  Semi-infinite  Medium 

B. l  General  Description 

The  SEMIl  program  (version  1)  has  been  developed  to  treat  neutral 
particle  transport  in  a  semiinfinite  medium.  The  analytical  solution  to  the 
transport  equation  for  single  and  two  half  space  geometries  for  a  variety  of 
sources  is  evaluated  numerically.  The  angular  flux  at  x  =  0  is  evaluated  us¬ 
ing  expressions  derived  from  Chandrasekhar’s  Principles  of  Invariance.  The 
scalar  flux  in  the  interior  is  evaluated  using  these  expressions  as  image  func¬ 
tions  for  a  Laplace  transform  inversion  performed  numerically  as  discussed 
above.  Six  source  configurations  axe  considered  for  each  half-space  geometry. 

B.2  Specific  Description 

B.2.1  Input 


The  input  to  the  program  is  as  follows; 


input 
line  1 

line  2 


description(free  format) 

number  of  cases 


line  3 


nc 

Ij 

Ig 

gt 


number  of  terms  in  inversion  series  (<100) 
number  of  iterations  of  bromwich  contour  (<10) 
to  determine  starting  contour  gamO  (0.1) 


nx 

iflux 


anl 

mx 


xO  initial  edit  position  from  x=0 

xl  final  edit  position  from  x=0 

number  of  intervals  between  xO  and  xl  {<m2) 

0  exiting  angular  flux  only 

1  interior  scalar  flux  only 

2  both  exiting  angular  and  scalar  fluxes 
limiting  direction  for  angular  edit  (-anl, anl) 
number  of  edit  points  (<ml) 

note:  repeat  line  4  for  the  nc  cases 

single  half-space  geometry: 
line  4  is  1  beam  source 

2  isotropic  source 

3  Milne  problem 

4  uniformly  distributed  isotropic  source 

5  distributed  isotropiexponential  source 

6  general  angular  distribution  at  surface 

two  half-spaces  geometry: 

51  beam  flux  at  interface 

52  half  range  source  at  interface 

53  Milne  problem 

54  uniformly  distributed  source 
in  right  half-space 

55  distributed  isotropic  exponential  source 
in  right  half-space 

56  general  angular  distribution  at  interface 

qO  source  normalization 

al  exponent  of  exponential  source 

be  exponent  for  general  angular  source 

(ff=e**(-be/(l-mu**2))) 
cl  number  of  secondaries  in  right-space 

c2  number  of  secondaries  in  left-space 

amO  source  direction 

Im  Gauss- Legendre  quadrature  order(<ml) 

err  desired  relative  error 

a)  Notes: 

1)  The  number  of  cases  nc  considered  is  unlimited. 
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2)  The  number  of  spatial  and  angular  edit  points  nx,inx  is  currently  lim¬ 
ited  to  100  for  the  single  half-space  and  200  for  two  half-spaces  but  can 
be  changed  by  specifying  ml  and  m2  in  the  program. 

3)  Ij  is  the  maximum  number  of  terms  allowed  in  the  inversion  series;  while 
Ig  specifies  the  maximum  number  of  contours  on  which  the  inversion  is 
evaluated  in  order  to  achieve  convergence. 

4)  gt  specifies  the  Bromwich  contour  7  =  7a  -hgt/a:  when  x  is  less  than  0.1 
and  is  the  initial  contour  otherwise.  7,  is  an  approximation  to  the  real 
part  of  the  rightmost  singularity. 

5)  The  general  angularly  dependent  source  currently  used  is 

/(//)  = 

where  ^  is  input  as  be  and  Qq  is  determined  from  normalization. 

6)  The  maximum  Gauss- Legendre  quadrature  order  is  ml 

B.2.2  Output 

The  exiting  or  interfacial  anguleir  flux  and  the  interior  scalar  flux  are 
output  on  files  tang  and  seal  respectively.  Diagnostic  messages  found  on  file 

o3  axe  as  follows: 

tape23  (o3)  diagnostic  messages: 

al  error  0:  h-function  iteration  did  not  converge 

b)  error  la,b:  search  for  zero  of  infinite  medium 

dispersion  relation  did  not  converge 

c)  error  2:  x=aaaaa  err=bbbbb: 

inversion  series  did  not  converge  at  x=aaaaa 
with  error  estimate  err=bbbbb 

d)  error  3:  possible  loss  of  accuracy  in  inversion  series 

e)  error  4:  j=iiiii  x=aaaaa  err=bbbbb: 

Romberg  integration  did  not  converge  at  term  j 
and  x=aaaaa  with  error  estimate  err=bbbbb 

f)  error  5:  x=aaaaa  err=bbbbb: 

contour  iteration  did  not  converge  at  x=aaaaa 
with  error  estimate  err=bbbbb 

g)  error  6:  search  for  zeros  of  the  Legendre  polynomial  failed 
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while  files  o4  and  o5  axe  plot  files  for  the  interior  scalar  flux  and  angular  flux 
at  X  =  0  respectively.  The  number  of  contours  for  convergence  is  monitored 
on  the  screen.  The  completion  of  each  case  is  signaled  by  an  indication  of 
the  files  which  have  been  written. 

B.2.3  Program  Notes 

1)  Loss  of  accturacy 

Since  the  inversion  is  performed  by  summing  a  series,  there  is  always 
the  possibility  of  loss  of  accuracy.  To  partially  offset  this  difficulty,  an 
iterative  procedure  has  been  introduced.  If  such  a  loss  is  sensed  by  the 
sum  being  less  than  any  one  term,  the  error  in  the  Romberg  integration 
is  reduced  to  obtain  more  accuracy  for  the  terms  of  the  series.  The 
computation  is  then  repeated. 

2)  Nonconvergence 

If  the  inversion  does  not  converge  for  small  r  (z  <  0.1),  gt  should  be  in¬ 
creased  in  an  attempt  to  achieve  convergence.  In  general,  the  inversion 
cannot  give  values  smaller  then  machine  accuracy.  Sometimes  by  judi¬ 
cious  a  choice  of  the  contour,  smaller  values  can  be  obtained  however. 
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MISSION 

OF 

ROMS  LABORATORY 


Roib«  Laboratory  plans  and  executes  an  interdisciplinary 
program  in  research,  development,  test,  and  technology 
transition  in  support  of  Air  Force  Command,  Control, 
Communications  and  Intelligence  (C3I)  activities  for  all 
Air  Force  platforms.  It  also  executes  selected 
acquisition  programs  in  several  areas  of  expertise. 
Technical  and  engineering  support  within  areas  of 
competence  is  provided  to  ESC  Program  Offices  <c\>s;  and 
other  ESC  elements  to  perform  effective  acquisition  of 
C3I  systems.  In  addition,  Rome  Laboratory's  technology 
supports  other  AFMC  Product  Divisions,  the  Air  Force  user 
community,  and  other  DOD  and  non-DOD  agencies.  Rome 
Laboratory  maintains  technical  competence  and  research 
programs  in  areas  including,  but  not  limited  to, 
communications,  command  and  control,  battle  management, 
intelligence  information  processing,  computational 
sciences  and  software  produclbility,  wide  area 
surveillance/sensors,  signal  processing,  solid  state 
sciences,  photonics,  electromagnetic  technology, 
superconductivity,  and  electronic 
reliability/maintainability  and  testability. 


